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Abstract. We search for evidence of dark matter (DM) annihilation in the isotropic 
gamma-ray background (IGRB) measured with 50 months of Fermi Large Area Tele¬ 
scope (LAT) observations. An improved theoretical description of the cosmological 
DM annihilation signal, based on two complementary techniques and assuming generic 
weakly interacting massive particle (WIMP) properties, renders more precise predic¬ 
tions compared to previous work. More specihcally, we estimate the cosmologically- 
induced gamma-ray intensity to have an uncertainty of a factor ~ 20 in canonical 
setups. We consistently include both the Galactic and extragalactic signals under the 
same theoretical framework, and study the impact of the former on the IGRB spectrum 
derivation. We hnd no evidence for a DM signal and we set limits on the DM-induced 
isotropic gamma-ray signal. Our limits are competitive for DM particle masses up to 
tens of TeV and, indeed, are the strongest limits derived from Fermi LAT data at TeV 
energies. This is possible thanks to the new Fermi LAT IGRB measurement, which 
now extends up to an energy of 820 GeV. We quantify uncertainties in detail and 
show the potential this type of search offers for testing the WIMP paradigm with a 
complementary and truly cosmological probe of DM particle signals. 

Keywords: gamma ray experiments, dark matter experiments, dark matter theory, 
gamma ray theory, dark matter simulations 
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1 Introduction 

The Fermi Large Area Telescope (LAT) [1] provides a unique potential to measure 
gamma-ray intensities with an almost uniform full-sky coverage in the energy range 
from 20 MeV to greater than 300 GeV. Due to its good angular resolution, more 
than 1800 gamma-ray point sources have been reported in the second source catalog 
(2FGL) [2] and more than 500 sources with a hard spectrum above 10 GeV have also 
been identihed in the high-energy IFHL catalog [3] . Most of these are of extragalactic 
origin. In addition, the excellent discrimination between charged particles and gamma 
rays allows LAT to directly measure diffuse gamma-ray emissions too. Note that this 
emission is notoriously hard to measure with Gherenkov telescopes from the ground at 
higher energies (above 100 GeV) due to isotropic cosmic-ray (GR) backgrounds (for a 
recent effort see [4]). As a result, the Fermi LAT is in a unique position to measure the 


1 




diffuse emission from the Milky Way with good angular resolution and to establish an 
isotropic emission, that is presumably of extragalactic origin, at energies greater than 
10 GeV [5], 

First detected by the SAS-2 satellite [6] and conhrmed by EGRET [7], the isotropic 
gamma-ray background (IGRB) is what remains of the extragalactic gamma-ray back¬ 
ground (EGB) after the contribution from the extragalactic sources detected in a given 
survey has been subtracted.^ The Fermi LAT collaboration has recently published a 
new measurement of the IGRB [8] based on 50 months of data and extending the 
analysis described in [5] down to 100 MeV and up to 820 GeV. The aim of this paper 
is to use this new measurement to search for evidence of a possible contribution from 
Weakly Interacting Massive Particle (WIMP) annihilation. This signal depends both 
on cosmological aspects of the DM clustering and the WIMP properties, and therefore 
potentially encodes a wealth of information. With the new measurements presented 
in ref. [8], it is possible to test DM models over a wide mass range, thereby testing 
candidates up to several tens of TeV (for Fermi LAT works that used other gamma-ray 
measurements for indirect DM searches see e.g., [9-11]). 

It is common to consider any DM annihilation signal viewed from the Sun’s po¬ 
sition as having three contributions with distinct morphological characteristics and 
spectra: Galactic smooth DM distribution. Galactic substructures and extragalactic 
(or, equivalently, cosmological) signal. The extragalactic DM annihilation signal is 
expected to be isotropic to a large degree and constitutes the main topic of our work, 
but we will carefully explore the relevance of the smooth Galactic and Galactic sub¬ 
structure components as well. This is one important addition to the methodology 
presented in the original Fermi LAT publication on this topic [12]. In addition to that 
we use improved theoretical predictions for the extragalactic DM signal, which both 
takes advantage of a better knowledge of the DM clustering and its cosmic evolution 
(section 2.1), and utilizes a complementary and novel method [13, 14] to calculate the 
extragalactic DM annihilation rate in Fourier space (section 2.2). Very interestingly, 
we hnd the complementary approach to agree well with the improved predictions from 
the traditional method (section 2.3). We also explore the degree of anisotropy of the 
DM signal that originates from Galactic substructures, and include that component 
consistently with the extragalactic DM emission (section 2.4). Finally, we present our 
constraints on the total isotropic DM signal in section 3, for which we make minimal 
assumptions on the isotropic signal of conventional astrophysical origin: we either as¬ 
sume this astrophysical emission to be negligible (and derive conservative limits) or 
we assume that its contribution is perfectly known and makes up the measured IGRB 
intensity (thereby estimating the sensitivity reach of the IGRB measurement to WIMP 
annihilation signals). We then study the robustness of the IGRB against adding a non- 

^ The EGB refers here to the high-latitude gamma-ray intensity after known diffuse Galactic 
contributions have been modeled and subtracted. Clearly, unmodeled Galactic components (including 
e.g. millisecond pulsars) would still be present in the EGB as well as in the IGRB. In particular, we 
will consider the Galactic DM halo with its subhalos as a possible contribution to the EGB — in spite 
of the fact that this might not be considered to be an “extragalactic” component. Moreover, emission 
from detected point sources is a known contributor and therefore we consider only the IGRB for our 
study of potential dark matter signals. 
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isotropic smooth Galactic DM halo signal in the IGRB derivation procednre of ref. [8]. 
Specihcally, we check for consistency of the IGRB measurement with the presence of 
such a Galactic DM signal in section 4. We summarize the main results of our work 
in section 5. 

2 Theoretical predictions of the isotropic dark matter anni¬ 
hilation signals 

The extragalactic gamma-ray intensity d(jP^/dEQ produced in annihilations of DM 
particles with mass and self-annihilation cross section (av) over all redshifts z is 
given by [15-17]: 

^0DM c{av){fl-[)MPc)‘^ f, + zYC{z) dN 

dEQ Svrm^ J H{z) dE 

where c is the speed of light, Ddm is the current DM abundance measured with respect 
to the critical density pc, H{z) is the Hubble parameter or expansion rate, and dN/dE 
is the spectrum of photons per DM annihilation.^ The function t{E, z) parametrizes 
the absorption of photons due to the extragalactic background light. The flux multiplier 
C{z), which is related to the variance of DM density in the Universe and measures the 
amount of DM clustering at each given redshift, is the most uncertain astrophysical 
quantity in this problem. It can be expressed both in real space, making use of a Halo 
Model (HM) approach [18], and in Fourier space by means of a Power Spectrum (PS) 
approach [13, 14]. 

In the HM framework, (C(z) is calculated by summing the contributions to the an¬ 
nihilation signal from individual halos^ of mass M from all cosmic redshifts, {E{M, z)), 
and for all halo masses, i.e.: 

where is the minimum halo mass considered, and A(z) and ^ are the DM 

halo-mass over-density and the halo mass function^ respectively. The mean halo over¬ 
density A{z) is defined with respect to the critical density of the Universe and its value 
determines the virial radius of a halo, Ra, at each redshift. In this paper we will adopt 
A(0) = 200, following previous choices in the literature, and compute it at different 
redshifts as in ref. [19]. The halo mass function ^ is normalized by imposing that 
all mass in the Universe resides inside halos (see [15] for more details). {E{M,z)) in 
turn depends on the DM halo density profile and the halo size. Halo density prohles 
are determined by N-body cosmological simulations, with the most recent results fa¬ 
voring cuspy Navarro-Frenk-White (NFW) [20] and Einasto halos [21, 22], while some 

^We assume here that the thermally-averaged annihilation cross section is velocity independent 
and that DM is composed of self-conjugated particles. 

^The term ‘halos’ refers to all types of virialized DM clumps and structures in the Universe that 
lead to a DM density enhancement over the background. 
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astrophysical observations favor cored halos, e.g., Burkert density profiles [23, 24], The 
density prohle k can be easily expressed in terms of a dimensionless variable x = r/vg, 
Tg being the scale radius at which the effective logarithmic slope of the prohle is —2. In 
this prescription, the virial radius is usually parametrized by the halo concentration 
ca = R/^/fg and the function F can be written as follows: 


F{M,z) 


c\{M,z) 


dx X^K^{x) 
[/q^ dx x"^ 


(2.3) 


More realistically F is an average over the probability distribution of the relevant pa¬ 
rameters (most notably ca). Note that the above expression depends on the third 
power of the concentration parameter. From simulations it has been determined that 
the halo mass function and halo concentration decrease with halo mass and, conse¬ 
quently, the hux multiplier (C(z) given by eq. (2.2) is dominated by small mass halos 
(which could be either held halos or subhalos^; see section 2.1 for further discussion). 
Furthermore, DM halos contain populations of subhalos, possibly characterized by dif¬ 
ferent mean values of the relevant parameters (e.g., diherent concentrations than those 
of held halos). The signals from subhalos are typically included by adding an extra 
term in eq. (2.2) to account for halo substructure, see [15]. 

As noted in [13] the hux multiplier can also be expressed directly in terms of 
the non-linear matter power spectrum P/vl (the Fourier transform of the two-point 
correlation function of the matter density held): 


c(z) ^ {nz)) = 


d k k^PNiik, z) _ 


k 


27r2 


^A^L{k,z), (2.4) 


where Aj^L{k,z) = k^P^i^k, z)/{27i^) is the dimensionless non-linear power spectrum 
and kma.x{z) is the scale of the smallest structures that signihcantly contribute to the 
cosmological annihilation signal. Loosely speaking, one could dehne a relation M = 
A/37rph {it/ kY with ph the characteristic density of the DM halo. Therefore fcmax is the 
PS counterpart to the minimal halo mass Mmin in eq. (2.2) in the HM prescription. 

The extrapolation to masses or k scales beyond the resolution of N-body simu¬ 
lations is the largest source of uncertainty in predictions of the extragalactic signal 
of DM annihilation, since the smallest scales expected for the WIMP models are far 
from being probed either by astrophysical observations or simulations.^ Thus, different 
methods of extrapolating to the smallest masses can lead to completely different re¬ 
sults for the relevant quantities. Typical expectations for the minimum halo masses in 
WIMP models are in the range Mmin G [10“®, lO”"^] Mq (see [28-30] and refs, therein), 
while we only have observational evidence of structures down to 10^ Mq [31] implying 

"^Subhalos are halos within the radius of another halo. A halo that does not reside inside any other 
halo will be referred to as a field halo or, simply, halo. 

^Notable exceptions on the simulation side are the works [25, 26], which simulated a few individual 
~ 1O“®M0 halos and the recent work [27], where for the first time dozens to thousands of halos 
have been resolved with superb mass resolution in the range 2 x 10“® — 4 x Note, though, 

that these simulations were performed at high redshifts, so some extrapolation is needed to describe 
structures at the present time. 
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that extrapolations are required to span > 10 orders of magnitude in halo mass (or 
> 3 orders of magnitude in k). 

Both ways of expressing (, eq. (2.2) and (2.4), have advantages and disadvantages. 
While eq. (2.2) is given in real space and thus deals with ‘intuitive’ quantities, it 
depends to a large extent on several poorly constrained parameters, most notably the 
concentration and halo mass function. This is particularly true for the smallest halos, 
which, as we have noted, are expected to dominate in the evaluation of (. The same 
is applicable to the subhalo population, whose internal properties and abundance are 
even less well understood. In contrast, eq. (2.4) depends only on one quantity directly 
determined from N-body simulations,® which can be extrapolated based on simple 
scale-invariant arguments, but lacks the intuitive interpretation of breaking structures 
into individual halos and subhalos, relevant when comparing the expected signals from 
Milky Way substructures with the total cosmological one. 

In this work, we use both of these approaches in parallel; the HM is used to dehne 
our benchmark model following simple but well motivated arguments for the choice of 
the relevant ingredients, and the PS framework is used to calculate the associated 
uncertainty due to extrapolation to small (unresolved) scales. 

2.1 Halo-model setup 

The DM annihilation signal from a halo depends on the third power of the concentration 
(see eq. 2.3), and the results will be extremely sensitive to the way the concentration- 
mass relation is extrapolated to low mass. A common practice in the past has been to 
use a single power law of M for c(M) all the way down to the minimum halo mass (see 
e.g. [32-34]). In most cases, these power-law extrapolations assign very high concentra¬ 
tions to the smallest halos, resulting in very high DM annihilation rates. In addition, 
these results are extremely sensitive to the power-law index used. However, these 
power-law extrapolations are unphysical and not expected in the ACDM cosmology 
[35]. Indeed, since natal CDM concentrations are set by the halo formation epoch and 
the smaller structures collapse at nearly the same time in the early Universe, low-mass 
halos are expected to possess similar natal concentrations, and therefore are expected 
to exhibit similar concentrations at the present epoch as well. In other words, the 
expectation that c(M) at the small mass end flattens is deeply rooted in the ACDM 
framework. This kind of behavior is correctly predicted by the c(M) model of ref. [36], 
which explicitly relates halo concentrations to the root mean square (r.m.s.) of the 
matter density fluctuations, a{M). In addition, both the ACDM expectations and the 
model predictions at the low-mass end are supported by the (few) results that come 
from N-body simulations that were specihcally designed to shed light on this extreme 
small-halo-mass regime [25-27]. This was also pointed out in ref. [35] where, making 
use of all available N-body simulation data, the authors examine the c(M) relation at 
redshift zero for all halo masses (i.e. from Earth-mass microhalos up to galaxy clus¬ 
ters). We refer the reader to these works for further details and discussion on the c(M) 
behavior at the low mass end. Note, though, that while the flattening of c(M) at lower 

®The quantity is the non-linear matter power spectrum, which is determined using only a matter 
density map, without invoking the concept of halos and without relying on standard halo finders. 
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masses is well motivated by ACDM and supported by simulations, the evolution of 
low-mass halos and subhalos from their formation times to the present and, in par¬ 
ticular, their survival probability, is far from being completely verihed in simulations 
or in analytical work (see, e.g.. Refs. [37-39]). Thus, even if the actual DM particle 
properties were known, unknowns in the subhalo survival probabilities and the impact 
of baryonic-feedback processes on the dark sector would still contribute to signihcant 
uncertainties on the theoretical expected DM signal. 

In what follows we dehne our benchmark halo model. 

• Cosmological parameters: we assume those recently derived from Planck data 
[40], i.e., Dm = 0.315, h = 0.673, ag = 0.834.^ 

• Minimum halo mass: here we make a choice that has become standard in the 
DM community, i.e. Mmin = 10“® h“^ Mq, though we caution the reader that in 
supersymmetric models this value is simply expected to be in the [10“®, 10“^] Mq 
range [28-30, 41]. We explore how predictions of the extragalactic DM signal are 
affected by the adopted Mmin value later in Section 2.3. 

• Halo concentration: we adopt the model of ref. [36] discussed above, which im¬ 
plicitly assumes a critical over-density, A in eq. (2.2), equal to 200; thus it gives 
a relation between C 200 and M 2 Q 0 . 

• Halo mass function: we use the state-of-the-art halo mass function as proposed 
by Tinker et al, [42], but with the parameters deduced by [36] at redshift zero 
for the Planck cosmology.® Inspired by the prior work by Sheth & Tormen [43], 
the Tinker et al. function gives a better agreement with N-body simulations 
especially at the high-mass end. 

• DM density prohle: we use the familiar NFW prohle [20]. The predicted values 
of the flux multiplier (C(z) are fairly independent of this choice (for the standard 
assumptions of the prohles considered in this work). It was shown for example, 
that the quantity {F{M, z)) in eq. (2.3) changes at the 10% level when the prohle 
is changed from the NFW prohle to the cored Burkert prohle (see also hgure 3 
in ref. [15]). 

• Contribution of the subhalo population: while both the halo mass function and 
halo concentrations are reasonably well studied, the properties of the subhalos 
are more uncertain. In order to estimate the DM annihilation rate produced in 
subhalos and its contribution to the total extragalactic signal, here we resort to 
the results of the study recently presented in ref. [35]. Two parameters that con¬ 
trol the amount of substructure and therefore its contribution to the annihilation 
hux are the minimum subhalo mass and the slope of the subhalo mass function. 
Following our choice for main DM halos, we adopt a value of 10“® Ii^^Mq for 

^The (Tg is the r.m.s. amplitude of linear matter fluctuations in 8 Mpc spheres at z = 0 with 
the Hubble constant h defined via H(z = 0) = h x 100 km/s/Mpc. 

®We still respect the parameters’ z-dependence found in ref. [42] though. 
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the first of these two ingredients. As for the slope of the subhalo mass function, 
oc {msub/Mhost)~°‘: we take a = 2, although we will also examine the impact 
of changing this exponent to 1.9 in the next section. Following results found in 
high-resolution N-body cosmological simulations of Milky-Way-sized halos above 
their mass resolution limits [44, 45] , the normalization of this subhalo mass func¬ 
tion is such that the mass contained in subhalos down to 10“® h~^MQ is ~ 45% 
of the total parent halo mass.® The mentioned values correspond to the refer¬ 
ence substructure boost model in ref. [35]. We note that this model implicitly 
assumes that both subhalos and held halos share the exact same internal prop¬ 
erties, which is probably not the case. Nevertheless, as discussed in ref. [35], this 
choice represents a conservative case in terms of expected gamma-ray intensity 
from annihilations in subhalos. 

We show the extragalactic intensity predicted by the benchmark HM described 
above in hgure 1. Since the uncertainty of this signal is not easy to quantify within the 
HM approach given its dependence on the many variables involved (see, e.g., [32, 46- 
48]), we use the PS approach to dehne our uncertainty band, as will be detailed in the 
next subsection. 

2.2 Power-spectrum setup 

In this section we estimate the uncertainty affecting the evaluation of the flux multi¬ 
plier, We focus on the PS approach, i.e. eq. (2.4), which relies on the work done in 
ref. [14] using the data from the Millenium-II N-body cosmological simulation (MS-II) 
[49]. MS-II has the highest mass resolution among available large-scale structure sim¬ 
ulations, but its cosmological parameters are from early WMAP 1-year [50] and 2dF 
Galaxy Redshift Survey [51] data. In ref. [14], it was shown that the largest uncertain¬ 
ties in the calculation of ( come from modeling the behavior of the power spectrum 
beyond the MS-II resolution and the exact value of the cut-off scale. Lacking a better 
theory, in [13, 14] the extrapolation of the data was inspired by the behavior of A^Lik) 
found in the simulation itself. In particular, at redshift zero, it was assumed that the 
true A 7 Vi(fc) value is bracketed by two alternative extrapolations for scales smaller 
than k > ki%\ 


A““(fc) = A^i(fci%) (2.5) 

A“^^(fc) = A^Lik,%) l—j (2.6) 

where Ueff = d In A(fci %)/dink and ki% is the scale at which the shot noise contribution 
to the power spectrum A^l is 1% and sets the resolution threshold. In other words, 
Amax(^) found by imposing that the spectral index Ues found at the resolution 
threshold stays constant at larger k scales. This is conservative, as a flattening of 

®Note, however, that the substructure mass fraction which is resolved in the simulations is only of 
about 10%. The 45% quoted in the text refers to the total (resolved plus unresolved) fraction which 
is obtained by extrapolating the subhalo mass function down to 10“® h^^M©. 


7 




the power spectrum is predicted by theoretical arguments (see, e.g., [14] for further 
details) while the extrapolation is from the trend observed in simulations at the smallest 
resolved mass scales. The clustering properties and the central parts of the DM density 
prohles at the smallest scales are not directly measured quantities and extrapolations 
are thus made under the assumptions mentioned above. 

At higher redshifts, where MS-II resolves only the largest co-moving scales, an 
additional constraint was imposed on the ratio of the nonlinear to linear 

power spectrum is not allowed to decrease with time at any hxed scale k (i.e., the 
Universe can only become clumpier in time, at every co-moving scale). Similarly, in 
order to also constrain the minimal ( value more precisely, A™™(A;) is required to have 
an effective spectral index bounded from below by the spectral index predicted in linear 
theory. For details on the physical motivations for these extrapolations see [14]. 

The cosmological parameters used for the MS-II simulation differ signihcantly 
from the values recently measured by Planck, which we adopted in the HM approach. 
The most critical factor in this regard is probably the value of ag, since this parameter 
has a large influence over the growth of fluctuations in the early Universe and thus on 
the subsequent evolution of structures. The results of the HM approach are derived 
assuming the most-up-to-date value ag = 0.835 recently given by Planck data, but 
MS-II adopts ag = 0.9. Thus, in order to make a fair comparison between the two 
methods, we follow [14] to apply the Planck cosmology to the results derived from 
the MS-II. We comment on the particular choice of the cut-off scale, fcmax; in the next 
subsection. 

2.3 Comparison of the two approaches and their cosmological dark matter 
signal predictions 

In hgure 1 we compare the values of ( obtained by means of the HM and PS approaches 
described above. In order to make a proper comparison and extract meaningful con¬ 
clusions, we call attention to a few caveats. 

A conceptual issue in comparing the results of the two approaches is the deter¬ 
mination of their cut-off scales. Such a cut-off exists, for example, due to the free 
streaming length of DM particles after kinetic decoupling [28, 29], which gives the 
highest frequency relevant Fourier mode fcmax-^° In the PS approach, fc max sets a sharp 
upper cut-off on the matter power spectrum at the scale corresponding to the presumed 
free streaming length in the linear regime, below which structures do not contribute 
to the DM annihilation signal. In the HM approach, the cut-off is instead imposed 
as a minimum halo mass,^^ below which no halos are formed and therefore no anni¬ 
hilation signal is expected. Within each DM halo, however, the annihilation signal is 
calculated by extrapolating the adopted DM density prohle of the halo down to r —)■ 0. 
Strictly speaking the HM thus includes Fourier modes all the way to inhnity, even if 
the largest wave numbers would not contribute much. Note that the exact shape of 
the DM density prohle at those small scales as well as the mass of the hrst virialized 
objects to form is still scarcely probed in simulations. As described in section 2.1 and 


the linear regime, the quantity /cmax is a redshift independent quantity. 
Chosen to correspond to the linear free streaming scale. 




z 


Figure 1. Normalized C as a function of redshift. A value of Mmin = 10“® Ii^^Mq was 
used in both the PS (gray) and HM predictions (red). The benchmark HM model detailed 
in section 2.1 is shown by the red solid line. The red dashed line corresponds to the case in 
which the slope of the subhalo mass function varies from the fiducial a = 2 to 1.9 (i.e., less 
substructure). The dotted line, labeled ‘PS (min)’, shows the minimum approximation from 
Equation 2.5 in the PS approach, while the dashed line, ‘PS (max)’, shows the maximum 
approximation given by Equation 2.6. 

2.2, we chose Mmin = 10“® h~^MQ as the cut-off scale in our benchmark HM model 
and we adopt fcmax = ^r/r^ as the default choice for the corresponding cut-off in the 
case of the PS approach (where is derived assuming the A = 200 value used in the 
HM approach, i.e. Vg = Ra{z)c/^{z) = R 2 oo(z)c 2 oo(z))- We recall that other motivated 
choices for fcmax are possible, and refer the reader to ref. [14] for further details. 

Another caveat when comparing HM to PS in figure 1 is the signal contribution 
from substructures that are present in extragalactic DM halos. As discussed in the 
previous section, the PS-derived ( values implicitly include such signal by covering 
contributions down to length scales ~ Tr/Zcmax- In the HM approach, by contrast, the 
substructures’ contribution is calculated by introducing additional parameters. We 
show in hgure 1 the HM prediction for two different scenarios: the ones correspond¬ 
ing to the minimum and maximum allowed values of the (substructure-induced) boost 
factor to the annihilation signal from held halos as predicted in ref. [35] for a hxed 
value of Mmin.ss = 10“® h“^M©. In this case, the differences in boost factors are due 
to different assumptions for the slope of the subhalo mass function, a (larger a val¬ 
ues lead to more substructure and thus to larger boosts). As a consequence of the 
aforementioned limitations, we expect some uncertainty when making a quantitative 
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Figure 2. Variation of the value of at z = 0 as a function of the minimum halo mass con¬ 
sidered, Mmin- Dotted and dashed lines represent minimum and maximum approximations 
in the PS approach as described by Eqs. (2.5) and (2.6), respectively. In this work, we adopt 
A^min = 10“® h“^ Mq as our fiducial value. See text for further details. 

comparison between the HM and PS approaches. Nevertheless, the agreement is qnite 
good as can be seen in hgnre 1, onr benchmark HM prediction being within the mini- 
mnm and maximum PS values at all redshifts. We have so far explored the expected 

WIMP signal for a given assumed cut-off scale Mmm (or, equivalently, fc max (^) dehned 
by tt/ts). However, this effective cut-off scale can vary signihcantly between various 
DM candidates, depending for example on their free-streaming lengths, as discussed 
in, e.g., [29]. In hgnre 2 we explore this dependence of ( on the cut-off scale and 

^max(^)- 

Finally, we now turn to the calculation of the extragalactic gamma-ray spectrum, 
by integrating ({z) over all redshifts and folding it with the induced spectrum from 
WIMP annihilations, see eq. (2.1). To model the attenuation of gamma rays traveling 
through cosmological distances, we adopt the Dominguez et al. model [52] for the 
extragalactic background light (EBL), which represents the state-of-the-art and is fully 
consistent with the hrst direct detections of the EBL attenuation by the LAT [53] and 
H.E.S.S. collaborations [54], and with the recently measured cosmic gamma-ray horizon 
[55]. In hgnre 3 we show a typical example for the gamma-ray hux resulting from 500 
GeV DM particles annihilating through the bb channel. As we will discuss in detail 
in section 3, the gray band in this hgnre translates directly into an uncertainty in the 
DM limits. 
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Figure 3. Comparison of the predicted cosmological DM-annihilation-induced gamma-ray 
intensities as given by both the PS and HM approaches. The C values implicitly used are 
those given in figure 1. This particular example is for a 500 GeV DM particle annihilating 
to bb channel with a cross section (av) = 3 x 10“^®cm^/s. The signal range shown in gray 
(which is computed within the PS approach) translates directly into uncertainties in the DM 
limits of section 3. 

2.4 Galactic dark matter signal contributions 

Within a DM halo, the DM distribution has two distinct components: a centrally- 
concentrated smooth distribution and a population of subhalos. While the smooth 
component has been studied extensively over a large span of halo masses, it is only more 
recently that several high-resolution simulations of Milky-Way-size halos addressed the 
subhalo component in more detail. These simulations quantified the radial distribution 
of subhalos inside the host halo, their abundance and their overall luminosity [44, 56, 
57]. The smooth and subhalo components are fundamentally different in terms of the 
level of anisotropy and, in the following, we quantify their relevance for this work. 

The smooth DM density profile of the main halo, p{r), is found in pure DM N-body 
simulations to be NFW [20] or Einasto [58, 59], while more computationally-expensive 
simulations which include baryonic effects have not yet converged on the profile shape in 
the inner regions of the Galaxy, finding both more cored [60, 61] and more cuspy profiles 
[62, 63]. Indeed, the solution to this issue might be substantially more complicated: 
the ability for galaxies to retain their DM cusps may depend on the ratio of their stellar 
and halo masses [64]. Observational tracers of the gravitational potential also cannot 
be used to determine the DM profile within < 2 kpc from the Galactic center as it is 
gravitationally dominated by baryons [65, 66]. However, the considerable uncertainties 
in the profile shape of the inner Galaxy are not critical for this work. We will deal 
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with Galactic latitudes ^ 20°, ~ 3 kpc from the Galactic Genter, a region in which 
simulations and astrophysical evidence converge on the p ~ profile behavior. 

Once the radial prohle of the DM density is hxed, the remaining uncertainty 
lies in its overall normalization, i.e. the value of po at the solar radius, which can 
take values in the range po = 0.2-0.8 GeVcm“^ [66, 67].^^ The gamma-ray signal is 
proportional to the square of p and therefore its uncertainty becomes greater than an 
order of magnitude in the worst case. We will not consider any (portion) of the signal 
from the smooth Galactic DM distribution to contribute to the isotropic emission (as 
has been done in some previous works, e.g. [68]). Instead, we find that the whole- 
sky DM template can be partially degenerate with at least one of the astrophysical 
components present in the Galactic foreground emission. We will treat this signal from 
the smooth Galactic DM halo as an additional component of the foreground Galactic 
diffuse emission instead. We further discuss this choice and its impact on the derivation 
of the IGRB spectra in section 4. 

As far as the diffuse gamma-ray intensity from DM annihilations in the Galactic 
subhalo population is concerned, some earlier works [12, 69] found that it could appear 
isotropic in our region of interest, and we explore this issue here in more detail. The 
exact distribution in Galactocentric distance of subhalos is currently not well deter¬ 
mined: in the Via Lactea II simulation [44] the subhalos follow the so-called anti-biased 
relation with respect to the smooth DM density profile, i.e. Psub(?") oc r x pNFw/(^" + ^"a) 
with Ta — 85.5 kpc, while in the Aquarius simulation [45, 56] they can be described 

with a particularly large scale radius 


by Psub oc PEinasto OC exp <{ - — 




^ - 1 


Tg = 199 kpc and ue = 0.678 [69]. 

In hgure 4 (left) we show that the substructures give on average < 10 % anisotropy 
in the relative intensity I / (/) for a DM-annihilation-induced signal when the DM 
substructure distribution is described as in the original Aquarius simulation paper 
[45]. However, for the same Aquarius simulation, authors in ref. [69] find a signal from 
the substructures that is less isotropic. 

Yet, the variations with respect to the average intensity are significantly less than 
a factor of 2 for latitudes \b\ > 20°. From the Via Lactea II simulation, the results 
in ref. [44] give an anisotropy of //(!) that is also less than a factor of 2, see hgure 4 
(right). These numbers can be compared to the Galactic smooth DM annihilation 
signal intensity that varies by more than a factor 30 for latitudes |6| > 20° for the 
NFW prohle that we used. 

We compare these hndings with the sky residuals found when deriving the IGRB 
spectrum in ref. [8], which are at the ^20% level. We conclude that, at least in the 
case of the Aquarius simulation, the Galactic substructure would lead to a sufficiently 
isotropic signal and thus we add it to the extragalactic signal when setting the DM 
limits. 

However, while the spatial distribution of the DM signal is taken directly from 


principle the values of and po are not independently constrained, and they should be corre¬ 
lated consistently. However, for our purposes the asserted ranges represent a fair description of the 
uncertainties expected for the DM signal. 
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Figure 4. Anisotropy in the gamma-ray annihilation signal from the subhalo distribution 
found in Aquarius [45] (left) and Via Lactea II [44] (right) simulations, with the former 
following the prescription in ref. [69]. The plots show the intensities of the substructures 
relative to their average intensity in the |6| > 20° region. 


the aforementioned simulations, the total signal strength is assigned self-consistently 
with the total substructure signal as of a typical Milky-Way-sized DM halo (as used 
in the calculation of the extragalactic signal within the HM approach). On the other 
hand, while the extragalactic DM signal comes from a large ensemble of halos, the 
properties of the Galactic substructure population are determined by the particular 
formation history of the Milky Way galaxy. In that sense, its properties do not have to 
correspond to the mean values found throughout the Universe, and in fact the Milky 
Way is found to be atypical in several respects [49]. For that reason, we consider two 
DM substructure prescriptions that introduce boosts of the total annihilation rate in 
our Galactic DM halo by factors of 3 and 15. This range follows from the prescription 
of [35] using a hxed minimum subhalo mass of Mmin,ss = 10“® h~^MQ but varying the 
slope of the subhalo mass function (a = 1.9 and 2, respectively). In the next section, we 
will show limits on DM annihilation cross sections from assuming these two bracketing 
values on the substructure boost. Ghanging «« to, e.g., 10“^^ h“^MQ would not 
affect the lower boost factor, but would increase the upper boost factor bound from 
15 to about 40 (see ref. [35]). 

Some of the largest or closest Galactic DM substructures could eventually be re¬ 
solved as discrete gamma-ray sources. The contribution from few individual subhalos to 
the total isotropic WIMP signal is not substantial, but nonetheless current constraints 
on DM signals from dwarf spheroidal galaxies [10, 11], as well as the non-detection of 
DM signals from unassociated gamma-ray sources, e.g., [70-72], signihcantly limit the 
total annihilation signal from the DM subhalos in the Milky Way. Our approach is 
to include the total expected DM signals from all subhalos of all masses in our eval¬ 
uation of the DM signal contribution to the IGRB, but when DM limits from, e.g., 
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dwarf spheroidal galaxies are stronger they obviously also impose limits on the total 
expected Galactic DM substructures contribution to the IGRB. Yet in these cases our 
limits are still relevant, as they represent an independent probe of cross sections by 
means of a conceptually different approach. 

While the gamma-ray signal originating from Galactic substructure could appear 
reasonably isotropic, an important difference with the extragalactic signal is in the 
spectral shape: the extragalactic signal is redshifted and distorted by absorption on 
the EBL (c.f. eq. (2.1)) while the Galactic signal directly reflects the injection spec¬ 
trum of gamma rays from DM annihilations and is generally harder. For that reason 
it is important to take this component properly into account, especially for heavy 
DM candidates, for which EBL absorption can severely limit high-energy gamma-ray 
intensities. 

3 Constraints on WIMP signals 

3.1 A review of the new Fermi LAT IGRB measurement 

Before we derive constraints on DM signals, let us summarize the four main steps that 
were taken in the analysis of ref. [8] to measure the IGRB, which will be used in this 
section and further discussed in the next section. In total, 50 months of LAT data 
were used, and dedicated cuts—creating two new event classes—were used to produce 
data samples with minimal contamination of GR backgrounds. In particular, the data 
were divided into a low-energy data set (events used to derive the spectrum in the 100 
MeV to 12.8 GeV range) and a high-energy sample (12.8 GeV to 820 GeV). Stronger 
cuts were applied in the low-energy range, where the event statistics are better but 
GR contamination is higher; the cuts were loosened at high energies, where both event 
statistics and GR contamination are lower. The full-sky data were then analyzed as 
follows: 

1. The full sky gamma-ray emission was modeled with a series of templates. Tem¬ 
plates of the Galactic diffuse emission are produced with the GALPROP code^^ 
[73], based on three distinct diffuse model setups, dubbed models A, B and G 
(see Appendix A and [8] for further details). The templates include maps of the 
gamma-ray emission due to interactions of GRs with interstellar gas and the in¬ 
verse Gompton (IG) emission separately. In addition, templates modeling the IG 
emission of GR electrons in the Solar radiation held, diffuse emission from Loop 
1 and point sources from the 2FGL catalog were used. After masking regions of 
bright interstellar emission along the Galactic plane, the normalization of each 
template was htted individually in energy bins in the range between 100 MeV 
and 12.8 GeV. 

2. Above 51.2 GeV the event statistics do not allow for htting in individual energy 
bins. To handle low statistics at high energies, the low-energy data set in the 
range from 6.4 GeV to 51.2 GeV was used to hnd the best-ht normalizations 

^^http: / / galprop.stanford.edu 
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of the Galactic templates. The normalizations were then hxed to these best- 
£t valnes, and GALPROP’s predictions of the spectral shapes were applied np to 
820 GeV in order to perform IGRB hts above 12.8 GeV using the high-energy 
data sample. 

3. The isotropic component thus derived is a sum of the IGRB emission and mis- 
classihed particle backgrounds. The IGRB is then obtained by subtracting a 
model for the GR contamination, obtained from Monte Garlo studies,from the 
isotropic emission. 

4. Using a baseline model for the Galactic diffuse emission (model A in ref. [8]), 
the above procedure was then repeated for different values of the relevant GR 
parameters. The scatter among the different IGRB spectra derived in this way, 
together with those derived by assuming foreground models B and G in [8], 
represents the systematic error band, indicated in hgure 12 of ref. [8]. 

The result of this procedure is a measure of the spectrum of the IGRB in the range 
from 100 MeV to 820 GeV [8]. It should be noted that this measurement is performed 
without including any Galactic smooth DM signal template, and the effects of such 
non-isotropic DM signal will be discussed in section 4. 

3.2 Known source contributors to the IGRB 

The Fermi LAT has detected many extragalactic sources: among the 1873 sources 
in the 2FGL catalog, there are 672 blazars (all classified according to the Roma BZ- 
GAT^^), 8 radio galaxies, 3 normal galaxies, 3 starburst galaxies and 2 Seyfert galaxies 
[2].^® The contribution to the IGRB from unresolved members of these extragalactic 
source classes has been studied over the years, e.g., [75-87]. In addition, some classes 
of Galactic sources, most notably millisecond pulsars, could contribute to the isotropic 
emission at large scale height in the Milky Way [88]. Their contribution to the IGRB 
is however severely constrained by the strong gamma-ray angular anisotropy signal 
expected for this source class [89, 90]. There are other truly diffuse emission processes 
that are expected to contribute to the IGRB as well, although probably only at a few 
percent level, e.g. structure formation shocks in clusters of galaxies [91] and giant radio 
lobes of FR II radio galaxies [92] . We refer the reader to Ref. [93] for a recent review 
on the IGRB and a discussion about potential contributors to this emission. 

Overall, the origin and composition of the IGRB are still open questions. Because 
of the large number of blazars detected by the LAT, direct population studies are now 
feasible using gamma rays and there is arguably a guaranteed contribution from the 

^^The Monte Carlo studies include a simulation of the relevant charged particle species and in¬ 
tensities present in the near-Earth environment as well as a phenomenological model for gamma-ray 
emission from the Earth limb. 

^®v4.1, http: // WWW .asdc.asi.it/bzcat/ 

^®There are 354 additional sources all associated in the 2FGL that appear to have blazar-like 
temporal or spectral characteristics but for which the lack of optical spectra did not allow a precise 
classification, most of them being labeled as AGN of uncertain type [74]. 
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blazar population [75-79]. The minimum contribution below 100 GeV from unresolved 
blazars has been estimated in ref. [76] to be close to 10%, the best estimate being 
22 — 34% in the 0.1-100 GeV range (which agrees well with previous hndings, e.g. 
[78, 94, 95]). The blazar contribution to the IGRB at the highest energies has only 
recently been studied. In ref. [78] they used a preliminary version of the new IGRB 
measurement reported in [8] and concluded that blazars can naturally explain the total 
measured IGRB above 100 GeV. 

For the other known source classes, however, we lack this kind of direct informa¬ 
tion, and cross correlations with radio (in the case of radio galaxies, see, e.g., [80-83]) 
or infrared data (for star-forming galaxies, e.g., [84]) have been used to determine the 
luminosity functions and infer the expected intensities in the Fermi LAT energy range. 
In a companion Fermi LAT paper [96], the contribution of blazars in the full energy 
range has been reevaluated using an updated luminosity function and spectral energy 
distribution model, taking advantage of recent follow-up observations [77] . When sum¬ 
ming the contribution from star-forming galaxies [84], radio galaxies [82] and blazars, 
ref. [96] shows that these three contributors could account for the intensity of the EGB 
across the 0.1 - 820 GeV range sampled by Fermi LAT. In ref. [96], the methodology 
of this work was adopted to derive DM limits taking advantage of the aforementioned 
new estimates of the astrophysical contributions. Yet, since intensity estimates for 
each of these potential IGRB contributors are uncertain (or under study) at the mo¬ 
ment, in this work we stay agnostic about the precise contribution of the astrophysical 
populations to the IGRB and instead aim for a more general approach. 

3.3 Statistical analysis 

We form a test statistic (TS) with a presumed distribution: 



(3.1) 


where Di is the measured IGRB intensity in energy bin i, V~^ is the inverse of the 
variance-covariance matrix and Mj is the IGRB model prediction (see below). 

On top of statistical uncertainties, the IGRB measurement inherits signihcant 
systematic uncertainties from the effective-area and the GR-contamination determina¬ 
tion. These systematic uncertainties, combined with the IGRB measurement procedure 
(summarized in section 3.1), can induce correlations between the IGRB measurements 
in the different energy bins. This should be reflected by a proper variance-covariance 
matrix V, and we made a study to estimate the variance-covariance matrix and check 
its impact on DM limits compared to the common approximation of taking it to be 
diagonal. 

To establish the expected correlation matrix for the data, 1000 Monte Garlo- 
generated pseudo experiments were created. Gamma-ray sky maps were generated 
with the help of HEALPix^' by taking the number of events in each sky pixel to be 
Poisson-distributed around the observed number of events. The effective area and 

^^http://healpix.jpl.nasa.gov/ 
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the CR contamination were drawn from normal distributions around their nominal 
values in [8]. To account for bin-to-bin correlations of the systematic uncertainty of 
the effective area we included correlations on the scale of three adjacent energy bins 
[97]. The CR contamination in the low- and high-energy data samples were taken 
to be fully uncorrelated. However, within each of these data samples, the systematic 
uncertainty (taken from [8]) was used to induce a randomized overall shift factor for the 
CR contamination rate level. The remaining (subdominant) statistical uncertainties 
for the CR contamination were taken to be fully uncorrelated. 

Each Monte Carlo-generated data set was then used to perform IGRB measure¬ 
ments exactly as done with the real data (using model A for the Galactic diffuse 
emission) and taking the GR and effective area determinations as described above. 
From the sample of 1000 IGRB pseudo measurements, the correlation matrix can then 
be directly calculated [98]. Subsequently, the variance-covariance matrix is determined 
by using the IGRB total variances (i.e. the sum of the variances from statistical, GR 
contamination and effective area uncertainties) in each energy bin as they were given 
in ref. [8]. The correlations are seen to be the strongest between neighboring energy 
bins at low energies, while energy bins further away and at the highest energies have 
negligible correlations. The derived variance-covariance matrix was then directly used 
in our statistical DM analysis (i.e. we included our calculated V in the TS calculation 
of eq. (3.1)), where it typically induced about a factor two increase in but less effects 
on Ay^. From this study we concluded that the impact on DM limits from including a 
proper variance-covariance matrix can have a sizeable effect but are typically smaller 
than the shift in DM limits coming from changing the diffuse foreground modeling 
shown in Appendix A. 

For our hnal analysis, we therefore decided to treat all data points D, as uncor¬ 
related with Gaussian probability distributions in our evaluations. The variance- 
covariance matrix V is thus approximated as diagonal with elements The system¬ 
atic uncertainties of the effective area and the charged GR contamination as well as 
the statistical errors are added in quadrature and their sum is the that enter in the 
covariance matrix V. The IGRB spectrum data points D* and the just specihed 1-cr 
errors cxj can be found in the supplementary tables of ref. [8].^® 

In addition, there is a signihcant systematic uncertainty due to the assumed Galac¬ 
tic foreground emission model in the IGRB derivation. The investigated range of 
Galactic foreground emission assumptions can induce correlated IGRB data points, 
and uncertainties are presumably very asymmetric. Galactic foreground uncertainties 
will therefore not be included in the evaluation of the y^ in eq. (3.1), but their impact 
was taken into account or estimated separately in our procedures, as we detail below. 

^®https://www-glast.stanford.edu/pub_data/845/ 
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3.3.1 WIMP signal search 

Before setting any limits, we will use a single power law with an exponential cut-off as 
a null-hypothesis background, 



(3.2) 


and search for any signihcant detection if a WIMP signal is added on top of this. This 
background should be viewed as the effective IGRB spectrum from all conventional as- 
trophysical sources discussed in section 3.2, where the overall normalization A, spectral 
index a and exponential cut-off energy Be are free parameters. 

For the DM signal search we use the IGRB derived with model A for the Galactic 
diffuse emission and only the Uj errors are included in the calculations. The simple 
background model of eq. (3.2) gives an excellent best-£t to the 26 data points of the 
IGRB spectrum; with a of 13.7 for 23 degrees-of-freedom.^® Naively this leaves little 
need for including an additional DM signal. To still search for a statistically significant 
DM signal we use the DM set-up described in section 2, where the isotropic DM signal 
is the sum of the contributions from cosmological signal and Galactic substructures. 
For the cosmological calculation we use the HM result, shown by the red solid line 
in hgure 1. Also, since the ratio of DM extragalactic to Galactic substructure signals 
affects the total DM signal spectral shape, we investigate the theoretical uncertainty 
range in the extragalactic signal strength by using the lowest and highest results from 
the PS approach (given by the limits of the gray band in figure 1). As for the Galactic 
substructures, in combination with either the HM or PS cosmological signals, we use 
two different overall signal strengths (corresponding to boosts of 3 and 15 for the Galaxy 
DM halo signal as a whole). The different WIMP annihilation channels we test are 
the same as those specihed in section 3.4, where we present our WIMP annihilation 
cross-section limits and derive our sensitivity reach. 

The x^ difference between the best-£t including such additional WIMP signal 
on top of the background (with the WIMP annihilation cross section as one extra 
degree of freedom) and the null hypothesis with zero DM signal, reveals that none of 
the DM hypotheses we tested showed a £t improvement by more than ATS = 8.3. 
Assuming that our TS has a x^ distribution, the local significance is 2.9a (before 
including any trials factor). The largest significance was in our minimal setup for both 
the extragalactic signal and the Galactic substructure signal, with a WIMP mass of 
500 GeV and an annihilation cross section into pairs of 1.1 x 10“^^ cm^/s. Note 
that this signihcance is not large, especially since a trials factor has not been yet applied 
(among 192 models we tested we had more than 2 ct detections for 16 different models). 
More importantly, uncertainties in the IGRB related to the selected Galactic diffuse 
emission model were not included. In fact, we also performed the analysis with IGRB 
derived with Galactic diffuse model B and G and conhrmed the DM non-detection 

^®This small ~xS value is presumably related to the fact that LAT’s systematic effective area and 
CR contamination uncertainties are included in the while their induced correlations are ignored in 
the variance-covariance matrix. If we include the variance-covariance matrix discussed in section 3.3 
then become 21.0; which is in good agreement with what should be expected. 
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obtained with model A. For example, the ATS = 8.3 mentioned above drops to 3.9 
when the IGRB derived with the Galactic diffuse model B is used. 

We test at most nine WIMP masses for each annihilation channel: 10, 50, 100, 500, 
1000, 5000, 10000, 20000, and 30000 GeV. This leaves a possibility that a somewhat 
more signihcant detection could be found if smaller steps in WIMP mass were used. 

For the tested WIMP masses we therefore conclude that there is no clear statis¬ 
tically signihcant evidence of WIMP signals in the IGRB, and we proceed to calculate 
upper limits. This is a rather naive approach and, indeed, a better understanding of 
the astrophysical contributions to the IGRB could help in revealing potential anomalies 
or even point toward the need for a DM contribution to the measured IGRB. 

3.3.2 Conservative approach for setting WIMP limits and the sensitivity 
reach of the IGRB measurement 

We will focus on deriving i) conservative DM limits, derived by making assumptions 
neither on the contributions from unresolved astrophysical populations to the IGRB 
nor on a specihc choice of a Galactic diffuse emission model, and ii) sensitivity reach, 
which assumes both that the total contribution from conventional astrophysical sources 
fully explains the measured IGRB at all energies by eq. (3.2), and that we can entirely 
rely on a specihc Galactic dihuse foreground model to derive the IGRB. 

We adopt the following two procedures: 

• Conservative limits: The 6Q- (3-1) is calculated only for bins where the 
DM signal alone exceeds the IGRB intensity: 




,2 


(3.3) 


cons 




where is the integrated DM-induced intensities in energy bin i as a function 
of (av) for a given WIMP candidate. Ehectively, this corresponds to having a 
(non-negative) background model that is free in normalization in each energy 
bin independently. To have a rough estimate of the ehect from having various 
Galactic foreground emission models that can alter the IGRB, we shift the centers 
of the IGRB data points derived with Galactic emission model A (while we keep 
the size of the errors cxj) to the upper edge of the 1-a envelope of all tested 
Galactic diffuse models used in ref. [8].^° This yields our data points The 

mentioned shift of the IGRB bins reflects the conservative approach of considering 
all tested diffuse emission models in ref. [8] equally likely for determining the 
IGRB (maximal) intensity in every energy bin. The 2a (3ct) DM limits are then 

We note that, in principle, one should marginalize over all the possible Galactic foreground 
models. Given that this is impractical, we opted here for a simpler approach that we believe provides 
a reasonable estimate of the effect. 
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defined to be the cases in which the DM signal component gives xLns eqnal to 4 


• Sensitivity reach: In this procednre, we use the IGRB derived with Galactic 
diffuse model A and include the DM signal on top of the background model in 
eq. (3.2). The of oq. (3.1) is then evaluated over all energy bins: 


Asens 


E 

i 




<T? 


(3.4) 


where (j)^^ + 0^'^® is the model prediction Mj as a function of {av) for a given 
WIMP candidate. Note that this represents a scenario in which i) the Galactic 
diffuse foreground used to derive the IGRB is hxed, ii) the contribution from 
conventional astrophysical sources to the IGRB is described by the parametric 
form of eq. (3.2) and hi) the parameters A, a, Ec in eq. (3.2) are fixed to their 
best-£t values found in ref. [8] (given in their table 4). The 2a (3(j) limits are 
then defined to be the cases in which the DM signal component forces the yfens 
to increase by more than 4 (9) with respect to the best-£t with a free DM 
signal normalization. 

We have also checked how this sensitivity reach changes by varying the adopted 
Galactic foreground model, namely by comparing limits when the IGRB is derived 
with the alternative foreground models A, B and G in ref. [8]. With this exercise, we 
gauge the impact of some systematic uncertainties associated with the modeling of the 
Galactic diffuse emission. We hnd differences that can be substantial especially for 
low WIMP masses; see Appendix A for further details. Yet, it should be noted that 
these tests are far from comprehensive and, as such, might not address the full range 
of uncertainties. 

The sensitivity reach derived here can also be taken as limits under the given 
assumptions. However, strictly speaking they should be interpreted as DM constraints 
only if the astrophysical background was independently predicted to the spectrum of 
eq. (3.2) with parameters equal to the best-£t values from the current IGRB measure¬ 
ment. 

The case where the total contribution to the IGRB from conventional astrophysics 
is derived as accurately as possible leads to DM constraints that typically lie between 
the conservative limit and the sensitivity reach derived in this work. Indeed, this is 
what is obtained in a companion work [96], where unresolved astrophysical source 
populations were modeled and used to set new DM limits on DM annihilation cross 
sections. 

avoid potential issues with correlations among the data points, we also performed the exercise 
of setting our limits by using one single bin (the one in which DM signal to the IGRB flux ratio is 
maximal) instead of a sum over bins as in eq. (3.3). We found that the two approaches lead to limits 
which differ by only E 10%. 
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Figure 5. Examples of DM-produced gamma-ray spectra which are at the border of being 
excluded by our 2a conservative limits. The WIMP mass and its annihilation channel is 
given in the upper left corner of each panel. The normalizations of the extragalactic signal 
and of the Galactic substructure signal are given by our benchmark HM model, as defined 
in section 2.1. Data points are in black, and the black lines show the upper and lower 
envelopes of the systematic uncertainties defined as the scatter among the different IGRB 
spectra derived in ref. [8]. 


In figures 5 and 6 we show illustrative examples of DM-induced spectra which 
have DM annihilation cross sections at the size of our 95% CL exclusion limits by our 
conservative approach and our sensitivity reaches, respectively. 

3.4 Limits on WIMP annihilation cross sections 

In this work, we stay agnostic about the nature of the DM particle and consider generic 
models in which DM annihilates with 100% branching ratio into bb, W~^W~, t^t~ or 
channels. For the hrst two channels, we consider only prompt emission and do 
not include any secondary gamma rays coming from the DM-induced electrons that 
up-scatter CMB photons by IC. Even for the heaviest DM masses we consider, the 
prompt emission is soft enough here to contribute signihcantly within the energy range 
of the measured IGRB - while the IC can be ignored because it only induces emission 
at much lower energies where the IGRB flux is higher. For and channels, 

instead, the prompt emission is harder and peaks signihcantly above the energy range 
for which the IGRB has been measured for our highest DM masses. In these cases the 
IC (which is also harder than for the previous two channels) contributes signihcantly 
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Figure 6. Examples of DM-produced gamma-ray spectra which are at the border of being 
excluded at 2a level in our procedure to calculate the sensitivity reach of the IGRB data. 
The WIMP mass and its annihilation channel is given in the upper left corner of each panel. 
The normalizations of the extragalactic signal and of the Galactic substructure signal are 
given by our benchmark HM model, as defined in section 2.1. Data points from ref. [8]. 


at energies close to the observed IGRB exponential cut-off and thus must be included. 
For that reason, both annihilation channels prove to be especially strongly constrained 
by the IGRB measurement [12]. We calculate the DM annihilation prompt spectra 
using the publicly available PPPC4DMID code [67], which takes into account electroweak 
bremsstrahlung corrections, which are particularly relevant for heavy DM candidates. 
For the calculation of the IG emission from the muon channel we follow the calculation 
presented in [12]. 

For the four annihilation channels under consideration, we present the conserva¬ 
tive limits and cross-section sensitivity reach at the 2a conhdence level in hgures 7 and 
8, respectively. In all cases, the DM limits were obtained by adopting the cosmologi¬ 
cal DM annihilation induced gamma-ray intensities given by the HM setup described 
in section 2.1, as well as a theoretical uncertainty range as estimated within the PS 
approach of section 2.2 (gray band in hgure 1). In addition, two conhgurations for the 
Galactic substructure contribution—which is assumed to be isotropic in this work—are 
adopted: i) the reference case, labeled as ’’SS-REF” in figures 7 and 8, where substruc¬ 
tures boost the total Galactic annihilation signal by a factor of 15, and ii) the minimal 
case, labeled ’’SS-MIN” in the figures, where the boost from Galactic substructure is 
3. Gonservative DM limits and cross-section sensitivities at the 3cr level for the bb and 
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T~^T channels were also derived, and can be found in Appendix B. 
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Figure 7. Upper limits (95% CL) on the DM annihilation cross section in our conservative 
procedure. From top to bottom and left to right, the limits are for the bb, W^W~, and 

channels. The red solid line shows limits obtained in our fiducial HM scenario described 
in section 2.1, and assumes the reference contribution from the Galactic subhalo population; 
see section 2.4 (‘HM, SS-REF’ case). The broad red band labeled as ‘PS (min—)-max), SS- 
REF’ shows the theoretical uncertainty in the extragalactic signal as given by the PS approach 
of section 2.2. The blue dashed line (‘HM, SS-MIN’) , with its corresponding uncertainty 
band (‘PS (min—i-max), SS-MIN’), refers instead to the limits obtained when the Milky 
Way substructure signal strength is taken to its lowest value as calculated in ref. [35]. For 
comparison, we also include other limits derived from observations with Fermi LAT [9, 11] 
and imaging air Cherenkov telescopes [99, 100]. 

From theoretical considerations, various DM particle candidate masses span a 
huge range. For thermally produced WIMPs, however, the Lee-Weinberg limit restricts 
the mass to be above few GeV [101] and unitarity considerations bound it to be below 
~ 100 TeV [102]. Interestingly, we are able to constrain signals for WIMP masses 
up to ~ 30 TeV because the IGRB measurement now extends up to 820 GeV. For 
DM particle masses above ~ 30 TeV, we start to probe the low-energy tail of the 
DM spectra and thus we lose constraining power rapidly. Furthermore, extragalactic 
WIMP signals are heavily suppressed at the highest energies as the optical depth is 
very large for such gamma rays. 

It is interesting to compare the conservative limits of hgure 7 to the cross-section 
sensitivities in hgure 8, at least for the case of our hducial HM scenario and the reference 
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Figure 8. DM annihilation cross section sensitivity reach (95% CL). Green solid line shows 
sensitivity obtained in our fiducial HM scenario described in section 2.1, and assumes the 
reference contribution from the Galactic subhalo population; see section 2.4 (‘HM, SS-REF’ 
case in the panels). The broad green band labeled as ‘PS (min—)-max), SS-REF’ shows 
the theoretical uncertainty in the extragalactic signal as given by the PS approach of sec¬ 
tion 2.2. The orange dashed line (‘HM, SS-MIN’), with its corresponding uncertainty band 
(‘PS (min—7-max), SS-MIN’), refers instead to the cross-section sensitivity obtained when the 
Milky Way substructure signal strength is taken to its lowest value as calculated in ref. [35] . 
For comparison, we also include other limits derived from observations with Fermi LAT [9, 11] 
and imaging air Cherenkov telescopes [99, 100]. 


contribution from the Galactic subhalo population (‘HM, SS-REF’ case in the panels). 
For the bb (r+r“) channel, the differences are of about factors 9, 25, 11, 3 (26, 9, 4, 3) 
at 10 GeV, 100 GeV, 1 TeV, 10 TeV. 

For low WIMP masses, the full spectral shape of the IGRB is affected by the 
WIMP signal, and hence the sensitivity reach, assuming a known spectral shape for 
the astrophysical contributions to the IGRB, places stronger limits, whereas for the 
largest WIMP masses only the last point (s) in the IGRB spectrum is affected and the 
two approaches are more similar. 

For the largest WIMP masses considered, the signal from Galactic substructures 

we omit the last data point, we find that both conservative limits and cross-section sensitivity 
for the bb channel worsen by < 30% at 5 TeV mass going up to a factor of --- 2 for masses between 10 
and 30 TeV. In the case of the harder channel, limits and sensitivity reach progressively weaken 
by a factor 2 to 4 between 2 and 30 TeV, respectively. 
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is stronger than that from the extragalactic DM, with the effect that the uncertainty 
range of the extragalactic WIMP signal becomes irrelevant when setting DM limits 
and calculating cross-section sensitivities. This is typically the case for gamma-ray 
energies above 100 GeV, where extragalactic signals are effectively attenuated due to 
EBL attenuation. The effect can be clearly seen in hgures 5 and 6 for several annihi¬ 
lation channels (see, e.g., the spectra of a 5 TeV mass DM particle annihilating to the 
W^W~). As a result, as the WIMP mass increases in hgures 7 and 8, the cross-section 
limit uncertainties get narrower (for a given Galactic substructure signal strength). For 
the same reason, the uncertainty band for the minimal Galactic substructure scenario 
(’SS-MIN’ case in hgures 7 and 8) is typically wider than the one for the reference 
Galactic substructure case (‘SS-REF’), especially at the largest WIMP masses consid¬ 
ered. This is less pronounced for the muon channel, because in that case the high-mass 
limits are still set by the IG peak of the emission which contributes at low energies. 

Another feature worth mentioning is that, in the case of DM annihilation into 
, hgures 7 and 8 show a dip in cross-section limits for DM particle masses around 
1 TeV. This dip is present because the part of the gamma-ray spectrum induced by 
FSR peaks at energies where the IGRB intensity has dropped exponentially above 
a few hundred GeV (see hgure 5). For larger WIMP masses, the FSR peak is well 
above the energy range covered by the LAT IGRB measurement and, as the WIMP 
mass increases, the limits get progressively weaker until lower-energy gamma rays— 
induced by IG upscattering of GMB photons from DM-induced high-energy electrons— 
eventually govern the constraints. 

Finally, it is also interesting to compare the conservative limits and cross-section 
sensitivities obtained in this work to other DM limits recently reported in the litera¬ 
ture. At low masses, ^ 100 GeV, the derived sensitivity reach is comparable to the 
limits derived from a stacking analysis of 25 dwarf spheroidal satellites of the Milky 
Way [11], and to the limits derived from considering diffuse emission at intermediate 
Galactic latitudes [9] . We note that the present analysis uses LAT data up to very high 
energies (820 GeV), which represents a novelty with respect to previous Fermi LAT 
DM searches. In order to put in perspective the results derived here, we also compare 
to DM constraints derived from ground-based atmospheric Gherenkov telescope obser¬ 
vations. Figures 7 and 8 show the limits derived by the Fl.E.S.S. Gollaboration from 
observing the Galactic center halo [99] and the MAGIC Gollaboration’s limits derived 
from deep observations of the Segue 1 dwarf spheroidal galaxy [100]. Our conservative 
limits are comparable to the latter ones in the TeV energy range, but weaker than those 
obtained by H.E.S.S.^^ As for the cross-section sensitivity reach, this is substantially 
better than the mentioned MAGIC limits, and comparable to the ones from H.E.S.S. 
Our work, which uses the IGRB’s total intensity to set constraints on the nature of 
DM, is connected to studies focused on small-scale angular anisotropies in the high- 
latitude gamma-ray sky [46, 103, 104]. Indeed, they both test for the presence of the 
same gamma-ray source class. In this regard the limits are directly comparable, and 

caution the reader that the H.E.S.S. limits are derived under the assumption of a cuspy profile, 
whereas DM limits from IGRB data are only moderately sensitive to the inner slope of the DM halo 
density profiles. 
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their DM signal limits are currently of similar magnitude [103, 104], It has also been 
shown that a cross-correlation of an anisotropic signal with the positions of galaxies 
in the 2MASS survey [105] or with weak lensing surveys [106, 107] could increase the 
ability to detect a DM component. 

4 Robustness of the IGRB measurement in the presence of a 
Galactic dark matter signal component 

The largest component of the systematic uncertainty in the measurement of the IGRB 
spectrum stems from the modeling of the Galactic diffuse emission [8]. A signal from 
the smooth Milky Way DM component would contribute to Galactic diffuse emission 
and, as seen in ref. [9], this signal can in part be degenerate with conventional as- 
trophysical emissions. Indeed, the DM component is morphologically similar to the 
inverse Gompton astrophysical emission and, in ref. [8], it is seen that uncertainties in 
the IG template have the most signihcant impact on the measured IGRB spectrum. 

In order to study this issue we repeat the htting procedure in [8] and model the 
Galactic diffuse emission, but with the addition of a Galactic smooth DM template. 
Our aim is twofold: i) to verify that DM and Galactic diffuse templates are partially 
degenerate with each other in the htting procedure, and ii) to a posteriori check for self- 
consistency of our procedure, i.e. for those DM annihilation cross sections constrained 
in section 3.4, we test whether the corresponding Galactic DM counterpart emission 
alters the IGRB measurement that was used to set the DM limits themselves. 

Following ref. [8], we perform low- and high-energy template hts separately (as 
described in section 3.1). We use templates of 10, 50, 250 and 500 GeV, and 1, 5, 
20 and 30 TeV DM particle masses annihilating to bb quarks and leptons, as 

representative gamma-ray signals from WIMP annihilation.^'^ We produce full-sky 
templates using the GALPROP code version 54 [108] into which we have incorporated 
the DM signal. We then evaluate the maximum value of the DM annihilation cross 
section (for each DM mass) which still leaves the IGRB measurement unchanged. The 
IGRB is taken as ‘unchanged’ when, after the inclusion of the Galactic smooth DM 
template, the new IGRB measurement for all energy bins falls inside either i) twice the 
width of the systematic uncertainty band derived from foreground model variations 
or ii) twice the la ‘statistical’ error bars of the IGRB measurement of ref. [8]. This 
maximal normalization can be translated into a maximal cross section once a particular 
Galactic DM density distribution is adopted. To be conservative—in the sense of 
hnding the corresponding DM annihilation cross-section values in our Galaxy which 
are ‘guaranteed’ to modify the IGRB—we set the local DM density to a low value of 
Po = 0.2 GeV cm“^ while keeping the Milky Way scale radius at the standard value of 
rg = 20 kpc, see e.g. [67].^^ 

^"*^Note that the morphology of the Galactic DM templates is independent of the DM annihilation 
channel considered, as long as one refers to the photons induced ‘directly’ in the annihilation process. 

^^Strictly speaking, rg and po Eire not independent, but for our purposes lowering po alone serves 
as a reasonable way to adopt a low but realistic Galactic DM signal. 
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Interestingly, while repeating the above procednre with different normalizations 
of the Galactic smooth DM templates we fonnd that, for DM masses below ~ 1 TeV, 
an increase of the Galactic smooth DM template normalization does not translate in a 
proportional reduction of the IGRB. Instead, this change gets compensated by a lower 
IG template which, as the DM normalization increases, can go down to almost zero 
in the corresponding energy bins before the IGRB measurement is altered. From this 
exercise, we conclude that a potential smooth Galactic DM signal could be (partially) 
modeled/absorbed by the IG templates, which thus may unintentionally subtract a DM 
signal in the IGRB measurement. This is also mostly the reason why we decided not to 
add any portion (e.g. the high-latitude emission) of the Galactic smooth DM template 
to the isotropic signal, as it was often done in literature, e.g. in refs. [68, 85]. For high- 
mass DM templates, however, degeneracies between the DM and the IG templates are 
limited, because the normalizations of the IG templates at these energies are hxed from 
gamma-ray data at intermediate energies (see point 2 in section 3.1). 

Note that this prescription does not guarantee that morphological residuals are 
kept small, since the methodology only investigates deviations in the IGRB spectrum. 
By also using the morphological information we could impose tighter constraints on a 
Galactic DM annihilation signal because large-scale anisotropic residuals could poten¬ 
tially start to appear in the IGRB measurement. In hgure 9, we show two examples of 
a changed IGRB spectrum for the two cases mentioned above in the case of a 1 TeV 
WIMP annihilating into bb. 

Figure 10 shows the largest possible DM annihilation cross sections to the bb and 
T^T~ channels which do not change the IGRB spectrum, together with our conservative 
limits on the cross section and sensitivity reach derived in section 3.^® The non-gray- 
shaded areas in hgure 10 roughly indicate the regions where our method of deriving 
limits on an isotropic DM signal would not lead to signihcantly altered results due to 
the modihed IGRB measurement from the presence of the assumed smooth Galactic 
DM signal. 

Notably, there are regions of the parameter space where DM limits overlap with 
the shaded areas of our conservative limits in hgure 10. Inclusion of the Galactic 
smooth DM template can lead to both smaller and larger IGRB intensities around 
the DM signal peak than the one reported in ref. [8]. For some DM masses ^ 250 
GeV the IGRB can e.g. get higher by up to ~ 40% after the inclusion of the DM 
template, which would naively weaken the limits by roughly this amount. For larger 
DM mass (^ 1 TeV) the IGRB spectrum is typically lowered. This is a consequence 
of our procedure in which the normalizations of the Galactic foreground spectra are 
determined at energies lower than the DM signal peak and then kept hxed at the higher 
energies. The measured IGRB intensity, the normalization of which is free in all energy 
bins, is therefore lowered around the energies where the DM signal peaks in order to 
accommodate the presence of the Galactic DM signal. 

Gross sections at the level of the sensitivity reach of the IGRB measurement are 
found to he below their gray shaded region. We check the impact of using modihed 

^®Only model A is used in the figure, but we note that similar results are obtained with models B 
and C. 


27 



10'^ 

^ .1 ^ 

^ .1 ^ ^ 


10"^ 

^ .1 ^ .1 ^ 

.1 ^ . 

7 

■ 

■ ■ 



7 

■ 

■ ■ 


cc 


■ ■ ■ . 


cc 



's 10-'' 


■ ■ 

1 

1 


'g 10-'' 





■ ■ 




■ I 



■ ■ ■ 


cc 


■ . ■ ■ 

> 


■ ■ 


> 



“ 10"'' 



— 

“ lO-'' 


■ ■ 

■ — 

-e- 




-e- 




. . original IGRB 




. . original IGRB 



■ ■ modified IGRB 




■ ■ modified IGRB 


10"® 

sys. error 

_,..1_^ 


10-® 

one sigma stat. error 

.1,, , , , 


10^ 10^ 10'' 10^ 10^ 10^ lo'* 10^ 


E [MeV] E [MeV] 

Figure 9. Left: The new IGRB measurement, after the inclusion of the Galactic smooth 
DM template, when the measurement for some energy bins falls outside twice the systematic 
uncertainty band, defined as the scatter among the different IGRB spectra derived in ref. [8] 
(the case to be compared with our conservative limits). Right: The modified IGRB, after 
the inclusion of the DM template, when the measurement for some energy bins falls outside 
two times the la statistical error band of the IGRB measurement originally presented in 
ref. [8] (to be compared with our calculation of the sensitivity reach). A 5 TeV WIMP 
which annihilates promptly into bb was used in both panels, which also explains why the 
maximum differences between the original and modified IGRB are found around 200 GeV in 
these particular examples. Note that we do not show the statistical error bar of the modified 
IGRB because it is not relevant for our determination of the modified IGRB. 

IGRB models derived under the assumptions that a Galactic DM signal is present. 
These alternate IGRB models are derived as above, with the Galactic DM signal fixed 
by the annihilation channel and cross section (the DM density profile is kept to the 
same as before). We adopt the cross-section values at the upper edge of the orange 
band in the top right panel of figure 10 (the ‘PS(min), SS-MIN’ case) and then apply 
our procedure to find the sensitivity reach: we find that the cross-section sensitivity 
curve is basically unchanged by the inclusion of the Galactic DM component. For cross 
sections within the gray shaded area the IGRB is sometimes no longer described well 
by the adopted background model, so the method is no longer expected to behave well. 

Note that while our DM limits depend on the substructure signal strength and 
the assumed minimal DM halo mass, the shaded gray region in hgure 10 is independent 
of it, so the relative position of the gray region and the limits would be different for a 
different choice of these parameters. 

In order to exhaustively explore the impact of Galactic smooth DM templates on 
the derivation of the IGRB, a larger number of Galactic astrophysical emission models 
should be studied. In this way it would be possible to probe in detail the IGRB along 
with various Galactic DM signals. However, such studies are beyond the scope of this 
work, which is tied to the methodology used in ref. [8]. The initial study performed in 
this section shows the importance of including the Galactic DM annihilation with its 
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Figure 10. The gray regions above the dotted lines indicate the DM annihilations cross 
sections which would alter the measured IGRB spectra significantly due to the signal from 
smooth DM halo component of the Milky Way; see section 4. Top and bottom panels are for 
bb and channels, respectively. The DM limits shown are the same as those presented 

in hgure 7 (left panels) and figure 8 (right panels). 


proper morphology in a detailed study of isotropic intensities.^^ Note that this issue 
is more severe for decaying DM models (studied in this context in e.g. [109, 110]), 
since the DM Galactic component is more isotropic compared to the annihilating DM 
case.^® Also, the Galactic DM signal at high latitudes is typically larger than the 
corresponding cosmological one [111], when compared to the annihilating DM case 
(see hgure 5), and therefore a careful study of the degeneracy between the DM and IG 
templates becomes mandatory before robust upper limits can be determined for any 
potential cosmological decaying DM signal. 

^^Note that the morphology of a DM template is relevant when performing the full-sky fits, since 
it influences the normalization of the isotropic template, i.e., the IGRB spectrum. This should not 
be confused, however, with the fact that we define the normalization of the DM template by the 
requirement that solely the IGRB spectrum gets changed, and not when the whole-sky residuals 
worsen significantly. 

^®More precisely, the smooth Galactic DM signal varies by factors of ^ 16 and ^ 4 for annihilation 
and decay, respectively, between Galactic latitudes 20°and 90°. 
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5 Summary 


We utilize the recent measurement of the IGRB spectrum, based on 50 months of 
Fermi LAT data, to set limits on the isotropic DM annihilation signals, i.e. the gamma 
rays originating from DM annihilation in halos over all of cosmic history as well as 
from the Galactic subhalos. Thanks to the broad energy range covered by the new 
IGRB measurement presented in ref. [8], which extends up to 820 GeV, we are able 
to effectively constrain signals from annihilation of DM particles with masses ranging 
between a few GeV and a few tens of TeV. 

At the lowest WIMP masses, our conservative DM limits in hgure 7 reach thermal 
cross-section values for bb and channels. For the case of our benchmark values for 
both the DM cosmological and Galactic substructure signals, the sensitivity reach of 
the IGRB measurement shown in hgure 8 is comparable to the limits recently obtained 
using LAT observations of dwarf spheroidal galaxies [11, 112] as well as those derived 
from low Galactic latitude data [9]. 

For WIMP masses above 5 TeV, our conservative limits calculated for the bench¬ 
mark values of the DM signal, are a factor of a few better than the ones presented in the 
Fermi LAT Gollaboration works cited above. At these high WIMP masses ( ^ 1 TeV), 
the benchmark conservative limits are comparable to those obtained from observations 
of dwarf spheroidal galaxies by ground-based gamma-ray telescopes (more precisely, 
the recent observation of Segue 1 by both VERITAS [113] and MAGIC [100]), but 
weaker than the limits derived from the Galactic center halo by H.E.S.S. [99]. The 
potential sensitivity to DM annihilation signals with the current IGRB measurement 
might reach an order of magnitude lower cross sections in this same WIMP mass 
range, in the case of optimal knowledge of some still uncertain (non DM) astrophysical 
factors. 

Our derived predictions for the strength of an isotropic DM annihilation signal 
can be considered realistic and not over estimated: the extrapolations performed below 
the resolution of current N-body cosmological simulations—necessary to account for 
the smallest halos—have been done in a physical and theoretically well-motivated way, 
and uncertainties in the expected DM signal were estimated using a well suited and 
complementary approach based on the non-linear matter power spectrum which is 
measured in N-body simulations. Furthermore, for the first time we have quantihed 
how the IGRB measurement is affected by a Galactic foreground DM signal and thus 
when the latter starts to impact the derived constraints on an isotropic DM signal. 

When compared to the earlier Fermi LAT Gollaboration work [12], which derived 
DM limits from the hrst-year IGRB measurement, the conservative limits are now 
about a factor of two stronger in the WIMP mass range 1 GeV to 1 TeV for the same 
value of the flux multiplier C(^)- This improvement can be attributed to the new 
IGRB data used and, most notably, to the fact that we did not take into account the 
(isotropic) signal from the Galactic substructure in the previous work. 

Moreover, the uncertainties of the flux multiplier ((z) have considerably shrunk 
in the present work, and it now has a factor ~ 20 uncertainty when the minimal 
halo mass cut-off is set to 10“® h“^ M©. We note that this theoretical uncertainty 
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range is a factor ~ 5 smaller than in ref. [12], We did not consider extreme power-law 
extrapolations of the many relevant qnantities in the Halo Model framework, which 
previously resulted in an over estimation of the predicted strength of the isotropic DM 
signal. In this work, the theoretical uncertainty range for the predicted DM signal 
strength (for a given WIMP annihilation cross section) therefore covers only the lower 
and more physically motivated part of the previously considered range in ref. [12]. 
This in turn implies that our limits are generally consistent with the most conservative 
estimates derived in other works [12, 85, 86]. 

In our work, we identihed and addressed three main sources of uncertainty affect¬ 
ing the derived limits to the DM annihilation cross section, which can be summarized 
as follows: i) theoretical predictions for the strength of the DM annihilation signal, 
which stem on one hand from the modeling of DM clustering at small scales, and which 
translate into an uncertainty of a factor ~ 20 (see section 2.3) and, on the other, from 
the precise amount of substructure in the Galaxy, which has a factor ~ 3 uncertainty, 
ii) modeling of the contribution of unresolved extragalactic sources to the IGRB, re¬ 
flected in the difference between our conservative limits and our derived sensitivity 
reach, which for bb and channels and our reference prediction of the isotropic 

DM annihilation signal ranges by a factor ~ 3 — 26, depending on the WIMP mass 
range considered, and iii) modeling of the Galactic diffuse emission, which can lead to 
variations in the limits by a factor of up to ~ 3. 

We also studied the impact on the IGRB measurement of a DM signal from the 
Milky Way DM halo, and dehned a region in the cross section versus DM mass plane 
for which the IGRB measurement would be sufficiently changed by the presence of 
Galactic DM as to potentially affect the limits derived here. 

With these considerable uncertainties in mind, at present the IGRB does not 
represent a clean target to search for a DM signal. At the same time, we showed in 
hgure 8 that the method has the great potential to test the ‘vanilla’ WIMP paradigm 
(i.e. the thermal cross-section value) up to masses of a few tens of GeV, making this 
approach competitive with other DM probes. An additional strength lies in offering 
a complementary and truly cosmological probe for any potential DM signal hint that 
might be claimed from another indirect, direct or collider search. 

In the coming years of the Fermi mission, the LAT sensitivity to point-like sources 
will improve due to the increased exposure and improved event classihcation and re¬ 
construction, which, depending on the energy band, will translate into a lower IGRB 
intensity and better DM limits. This can be the case at energies above 100 GeV, 
where the total IGRB might ultimately be attributable to Fermi LAT-detected point¬ 
like sources (as suggested in ref. [8], where the cumulative intensity from 2FGL sources 
and the IGRB are seen to be comparable at the highest energies). On the other hand, 
at energies below 100 GeV, resolving more point sources in the next several years is not 
expected to signihcantly lower the IGRB flux (since, e.g., new blazars will be extremely 
faint and thus their contribution will be hard to detect by the LAT, see [76]). Most 
importantly, developments are foreseen to help in lowering the most critical uncertain¬ 
ties: i) the Euclid satellite [114] and next generation N-body cosmological simulations 
(e.g. [115, 116]), by shedding light on the small-scale DM clustering properties, ii) new 
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constraints on the contribntion of astrophysical sources to IGRB, from anisotropy and 
cross correlation studies of high latitude gamma-ray emission, and iii) more precise 
cosmic-ray measurements (with e.g. AMS-02 [117]) as well as detailed Planck dust 
maps, that will help in rehning the modeling of the Galactic foreground. 
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Some of the results in this paper have been derived using the HEALPix code [118]. 

A Diffuse foreground models and their impact on limits 

In this Appendix, we investigate how the limits on the DM annihilation cross section 
depend on the particular foreground model that is used to derive the IGRB in ref. [8]. 
In section 3.4 the baseline foreground emission model A in ref. [8] was assumed for 
computing the DM limits. However, two more foreground models, B and G, were also 
dehned and considered in [8], which led to slightly different derived spectra of the 
IGRB. These three reference models differ in propagation and GR injection scenarios, 
and have been derived from a customized version of GALPROP. Model A is the basic 
reference model, whereas model B includes, e.g., an additional population of electron- 
only sources near the Galactic center and model G allows the GR diffusion and re¬ 
acceleration to vary signihcantly throughout the Galaxy. Furthermore, variations of 
model A have been studied, by, e.g., changing the size of GR halo between 4 kpc and 
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10 kpc, modifying its CR-source distribution, turning off re-acceleration and adding a 
‘Fermi Bubbles’ template to assess the systematic uncertainties in the derived IGRB. 

In hgure 11, we show the effect on the limits when assuming models B and C 
for the foreground diffuse emission instead of model A, for the particular case of our 
sensitivity reach estimate and the bb and r+r” annihilation channels. 

The limits are substantially modihed for low mass WIMPs, i.e., below ~ 300 GeV, 
although the limits are again rather comparable at the lowest DM masses considered, 
namely below ~ 20 GeV. The maximum difference between the limits is found at about 
50 GeV, where model A makes the limits a factor ~ 2.7 more stringent than the ones 
deduced using model B (assuming our hducial Galactic substructure scenario). In all 
cases, these differences can be explained as the interplay between the measured IGRB 
points obtained in ref. [8] for each foreground model and the WIMP mass considered. 
For example, in the particular case of a 50 GeV WIMP, for which the emission peaks 
at roughly a few GeV, there is a downward trend of the IGRB data between 2 GeV 
and 10 GeV in models A and G, which is not present in model B, making the limits 
substantially stronger for models A and G compared to B. 

It is worth emphasizing that, in contrast, the conservative limits ‘by construction’ 
take into account the variation induced on the IGRB measurement from using the 
different foreground models. In this case, we recall that, in order to set the limits, we 
shifted the IGRB data points to the maximum allowed intensity values among those 
given by the various foreground models, which always results in the most conservative 
bounds. 
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Figure 11. DM limits for the three IGRB measurements derived assuming foreground 
emission models A, B and C in ref. [8]. The limits shown in section 3.4 implicitly assume 
foreground model A in ref. [8] . The limits are substantially modified when assuming model B 
instead, especially at low masses. This figure is for the particular case of the sensitivity reach 
procedure (as it was described in section 3.3.2) for DM particles annihilating into bb quarks 
(top) and T~^T~ {bottom), and for the two scenarios of the Galactic substructure contribution 
introduced in section 2.4. The case of our reference substructure model is shown on the left 
panels, while the minimal substructure case is on the right. 


34 
























B Limits at different confidence levels 


In this Appendix we compare 2a and 3cr upper limits on the DM annihilation cross 
section for bb and channels (hgures 12 and 13, respectively), for both the conser¬ 
vative limits and the sensitivity reach, and for the six representative cases of the DM 
signal strength considered in our work. 
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Figure 12. Upper limits on the self-annihilation cross section obtained in our conservative 
(left) and sensitivity reach (right) procedure. The limits are on the annihilation cross section 
into bb quarks and the thicker (thinner) lines show the 2cj (Su) limits. The figures in each 
row show the limits in different DM setups: (top) the minimal extragalactic signal in the PS 
approach, (middle) the benchmark extragalactic signal in the HM approach, and (bottom) the 
maximal extragalactic signal in the PS approach. In each figure, the solid line corresponds to 
the benchmark Galactic substructure intensity, while the dashed line represents the minimal 
assumed Galactic substructure signal; see section 2.4. The minimal scale for DM structures 
corresponding always to a halo mass cut-off of Mmin = 10“® h“^M 0 . 
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Figure 13. Same as figure 12, but for the r+r 


channel. 
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